# JCR REPLICATION for "How Civilian Loyalties Shape Rebel-led Victimization of Rebel Constituencies" #
# ILAYDA B. ONDER #
# DESCRIPTIVE STATISTICS, DESCRIPTIVE FIGURES (MAPS), BALANCE TESTS #
# MAY 27, 2024 #
# This .R  file generates Table 2, Figure 1, Figure 2, Table 3, and Figure 4

#install.packages("dplyr")
library(dplyr)

# TABLE 2----
load("Descriptives Data.RData") # load descriptive statistics data

# Restrict the dataset to the variables used in the main analysis
df_sum = df_descriptives[c('publicworker', 'collaborator', 'politician', 'margin_kurd', 'govrep_vil', 'govrep_non', 'frequencyclash', 'intensityclash', 'murders', 'deaths', 'voters', 'recruitswe', 
               'pop', 'urbanpop', 'kurdvote', 'kurdpop', 'distance', 'altitude', 'rebellion', 'literacy', 'soceco')]
colnames(df_sum)[4] <- "marginkurd"
colnames(df_sum)[5] <- "govrepvil"
colnames(df_sum)[6] <- "govrepnon"

df_sum2 <- df_sum %>% # Summary by group using dplyr
  summarise_all(funs(mean, min, max, sd), na.rm = TRUE) %>%
  tidyr::pivot_longer(cols = everything(), names_sep = "_", names_to  = c("Variable", ".value"))
colnames(df_sum2)[2] <- "Mean"
colnames(df_sum2)[3] <- "Min"
colnames(df_sum2)[4] <- "Max"
colnames(df_sum2)[5] <- "SD"
df_sum2[2:5] <- round(df_sum2[2:5], 2)

# Rename columns according to variable names
df_sum2[1] = c("Public Workers (number of incidents)", 
               "Traitors (number of incidents)",
               "Pro-Government Local Politicians (number of incidents)",
               "Incumbent Party Vote Margin",
               "Violent Government Repression (number of incidents)",
               "Non-Violent Government Repression (number of incidents)",
               "Frequency of armed clashes prior to 2012", 
               "Casualties from armed clashes prior to 2012", 
               "Extrajudicial killings/political assassinations prior to 2004", 
               "Insurgent casualties prior to 2012", 
               "Number of voters in 2014 Municipal Elections", 
               "Insurgent recruits prior to 2012", 
               "District population", 
               "Urban population", 
               "Kurdish political party vote share in elections prior to 2014", 
               "% of Kurdish population", 
               "Distance to capital city", 
               "Altitude of district centers", 
               "Rebellion in the 1920s", 
               "Literacy rate", 
               "Level of socioeconomic development")
print(df_sum2, n = 21) # Table 2
View(df_sum2) # check df_sum2 with the command "View(df_sum2)" to see the exact numbers with decimal points





# FIGURE 1 ----
#install.packages("tidyverse")
#install.packages("sp")
#install.packages("ggplot2")

library(tidyverse)
library(sp) # package for spatial data
library(ggplot2)

load("Map Data DVs.RData") # load map data for DVs

map1 = ggplot(df_mapDV) +
  geom_polygon(aes(x = long, y = lat, group = group, fill = publicworker), color = "black") +
  coord_map() +
  theme_void() + 
  binned_scale(aesthetics = "fill",
               scale_name = "stepsn", 
               palette = function(x) c("seashell1", "pink", "hotpink2", "firebrick2", "firebrick"),
               breaks = c(1, 2, 3, 10),
               limits = c(0, 10),
               show.limits = TRUE, 
               guide = "colorsteps",
               na.value = "#8c8c8c") +
  theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.height =  unit(0.5, 'cm'),
        legend.key.width =  unit(4, 'cm'),
        legend.text = element_text(size = 18)) # Generate Panel A of Figure 1
map1

map2 = ggplot(df_mapDV) +
  geom_polygon(aes(x = long, y = lat, group = group, fill = collaborator), color = "black") +
  coord_map() +
  theme_void() + 
  binned_scale(aesthetics = "fill",
               scale_name = "stepsn", 
               palette = function(x) c("seashell1", "pink", "hotpink2", "firebrick2", "firebrick"),
               breaks = c(1, 2, 3, 10),
               limits = c(0, 10),
               show.limits = TRUE, 
               guide = "colorsteps",
               na.value = "#8c8c8c") +
  theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.height =  unit(0.5, 'cm'),
        legend.key.width =  unit(4, 'cm'),
        legend.text = element_text(size = 18)) # Generate Panel B of Figure 1
map2

map3 = ggplot(df_mapDV) +
  geom_polygon(aes(x = long, y = lat, group = group, fill = politician), color = "black") +
  coord_map() +
  theme_void() + 
  binned_scale(aesthetics = "fill",
               scale_name = "stepsn", 
               palette = function(x) c("seashell1", "pink", "hotpink2", "firebrick2", "firebrick"),
               breaks = c(1, 2, 3, 10),
               limits = c(0, 10),
               show.limits = TRUE, 
               guide = "colorsteps",
               na.value = "#8c8c8c") +
  theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.height =  unit(0.5, 'cm'),
        legend.key.width =  unit(4, 'cm'),
        legend.text = element_text(size = 18)) # Generate Panel C of Figure 1
map3






# FIGURE 2----
#install.packages("tidyverse")
#install.packages("sp")
#install.packages("ggplot2")

library(tidyverse)
library(sp) # package for spatial data
library(ggplot2)

load("Map Data IVs.RData") # load map data for IVs

map4 = ggplot(df_mapIV) +
  geom_polygon(aes(x = long, y = lat, group = group, fill = margin_kurd), color = "black") +
  coord_map() +
  theme_void() + 
  scale_fill_gradient(low = "#DC0000B2", high = "#4DBBD5B2") +
  theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.height =  unit(0.5, 'cm'),
        legend.key.width =  unit(4, 'cm'),
        legend.text = element_text(size = 18)) # Generate Panel A of Figure 2
map4

map5 = ggplot(df_mapIV) +
  geom_polygon(aes(x = long, y = lat, group = group, fill = akp_win), color = "black") +
  coord_map() +
  theme_void() + 
  scale_fill_manual(values=c("#DC0000B2", "#4DBBD5B2", "#8c8c8c")) + 
  theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.height =  unit(0.5, 'cm'),
        legend.key.width =  unit(4, 'cm'),
        legend.text = element_text(size = 18)) # Generate Panel B of Figure 2
map5

map6 = ggplot(df_mapIV) +
  geom_polygon(aes(x = long, y = lat, group = group, fill = close_akp_win), color = "black") +
  coord_map() +
  theme_void() + 
  scale_fill_manual(values=c("#DC0000B2", "#4DBBD5B2", "#8c8c8c", "seashell1")) + 
  theme(plot.margin = unit(c(0.2,0.2,0.9,0.2), "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.height =  unit(0.5, 'cm'),
        legend.key.width =  unit(4, 'cm'),
        legend.text = element_text(size = 16)) # Generate Panel C of Figure 2
map6





# TABLE 3----
load("Balance Test Data.RData") # load data for balance tests
sum = df_balance %>%  # compute means of variables
  group_by(akp_win) %>% # group by treatment variable
  summarize(govrep_vil = mean(govrep_vil, na.rm = TRUE),
            govrep_non = mean(govrep_non, na.rm = T),
            frequencyclash = mean(frequencyclash, na.rm = TRUE), # Frequency of armed clashes prior to 2012
            intensityclash = mean(intensityclash, na.rm = TRUE), # Casualties from armed clashes prior to 2012
            murders = mean(murders, na.rm = TRUE), # Extrajudicial killings/political assassinations prior to 2004
            deaths = mean(deaths, na.rm = TRUE), # Insurgent Casualties
            recruitswe = mean(recruitswe, na.rm = TRUE), # Insurgent recruits
            voters = mean(voters, na.rm = TRUE), # Number of voters in 2014 Municipal Elections
            pop = mean(pop, na.rm = TRUE), # District population
            urbanpop = mean(urbanpop, na.rm = TRUE), # Urban population
            kurdvote = mean(kurdvote, na.rm = TRUE), # Kurdish political party voteshare in elections prior to 2014
            kurdpop = mean(kurdpop, na.rm = TRUE), # % of Kurdish population
            distance = mean(distance, na.rm = TRUE), # Distance to capital city
            altitude = mean(altitude, na.rm = TRUE), # Altitude of district centers
            rebellion = mean(rebellion, na.rm = TRUE), # Rebellion in the 1920s
            literacy = mean(literacy, na.rm = TRUE), # Literacy rate
            soceco = mean(soceco, na.rm = TRUE)) # Level of socioeconomic development # International border

# Run t-tests
ttests <- lapply(df_balance[,c('govrep_vil', 'govrep_non',
                        'frequencyclash', 'intensityclash', 'murders', 'deaths', 'recruitswe', 'voters', 
                        'pop', 'urbanpop', 'kurdvote', 'kurdpop', 'distance', 'altitude', 'rebellion', 'literacy', 'soceco')], 
                 function(x) t.test(x ~ df_balance$akp_win, var.equal = TRUE))


# rearrange the table
sum2 = data.frame(t(sum))
colnames(sum2)[1] <- "Incumbent Party Loss"
colnames(sum2)[2] <- "Incumbent Party Victory"
sum2 <- sum2[-1, ]
sum2$ttest = c(ttests$govrep_vil$p.value, # add t-test statistics to the table
               ttests$govrep_non$p.value,
               ttests$frequencyclash$p.value,
               ttests$intensityclash$p.value,
               ttests$murders$p.value,
               ttests$deaths$p.value,
               ttests$recruitswe$p.value,
               ttests$voters$p.value,
               ttests$pop$p.value,
               ttests$urbanpop$p.value,
               ttests$kurdvote$p.value,
               ttests$kurdpop$p.value,
               ttests$distance$p.value,
               ttests$altitude$p.value,
               ttests$rebellion$p.value,
               ttests$literacy$p.value,
               ttests$soceco$p.value)
sum2$ttest = round(sum2$ttest, 2)

# set row names
rownames(sum2) <- c("Violent Government Repression (number of incidents)",
                    "Non-Violent Government Repression (number of incidents)",
                    "Frequency of armed clashes prior to 2012", 
                    "Casualties from armed clashes prior to 2012", 
                    "Extrajudicial killings/political assassinations prior to 2004", 
                    "Insurgent casualties prior to 2012", 
                    "Insurgent recruits prior to 2012", 
                    "Number of voters in 2014 Municipal Elections", 
                    "District population", 
                    "Urban population", 
                    "Kurdish political party vote share in elections prior to 2014", 
                    "% of Kurdish population", 
                    "Distance to capital city", 
                    "Altitude of district centers", 
                    "Rebellion in the 1920s", 
                    "Literacy rate", 
                    "Level of socioeconomic development")
sum2$`Incumbent Party Loss` = round(as.numeric(sum2$`Incumbent Party Loss`), 2) # rounding
sum2$`Incumbent Party Victory`= round(as.numeric(sum2$`Incumbent Party Victory`), 2) # rounding
sum2 # generate Table 3




# FIGURE 4----
load("Descriptives Data.RData") # Load descriptives data

#install.packages("rdd")
library(rdd)
options(scipen = 999)

par(cex.axis = 1, cex.lab = 1, cex.main = 1, cex.sub = 1) # figure configurations
par(mgp = c(3,1,0), mar = c(2, 2, 2, 2)) # figure configurations
par(mfrow = c(4, 4)) # figure configurations

plot(RDestimate(govrep_vil ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Violent Government Repression")

plot(RDestimate(govrep_non ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Non-Violent Government Repression")

plot(RDestimate(frequencyclash ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Armed Clashes")

plot(RDestimate(intensityclash ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Casualties from Armed Clashes")

plot(RDestimate(murders ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Extrajudicial Killings")

plot(RDestimate(deaths ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Insurgent Casualties")

plot(RDestimate(recruitswe ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Insurgent Recruits")

plot(RDestimate(voters ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Number of Voters")

plot(RDestimate(pop ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "District Population")

plot(RDestimate(urbanpop ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Urban Population")

plot(RDestimate(kurdvote ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Kurdish Party Voteshare prior to 2014")

plot(RDestimate(kurdpop ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "% of Kurdish Population")

plot(RDestimate(distance ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Distance to Capital")

plot(RDestimate(altitude ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Altitude")

plot(RDestimate(rebellion ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Rebellion in the 1920s")

plot(RDestimate(literacy ~ margin_kurd, data = df_descriptives, cutpoint = 0))
abline(v = 0, col = "red")
title(xlab = 'Incumbent Party Vote Margin', 
      main = "Literacy Rate")






